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SECTION  I 


INTRODUCTION  AND  SUMMARY 

This  special  report  investigates  several  practical  aspects  of 
generating  high- resolution  wavenumber  spectra  using  subarray  outputs  of 
the  Montana  LASA.  The  following  are  questions  of  particular  interest: 

•  Variability  of  traveltime  anomalies  as  a 
function  of  wavenumber 

•  Spectral  window  effect  on  crosspower  esti¬ 
mates  due  to  moveout  across  the  array 

•  Tradeoffs  involved  in  a  finite  length  trans¬ 
form  of  array  data 

These  questions  are  investigated  with  either  easily  available 
data  or  simple  models  of  the  situation  being  considered.  A  general  under¬ 
standing  of  the  type  of  constraints  encountered,  rather  than  an  exact  mathe¬ 
matical  treatment  of  each  situation,  is  desired. 

A.  TRAVELTIME  ANOMALY  VARIATIONS 

Variations  in  traveltime  anomalies  were  analyzed  using 
punched  cards  obtained  from  the  Seismic  Data  Laboratory,  Alexandria,  Vir¬ 
ginia.  The  punched  cards  contained  the  date,  time,  latitude,  longitude, 
and  arrival  time  recorded  at  each  LASA  subarray  for  384  teleseismic  events. 

A  map  of  the  world,  as  seen  teleseismically  (in  k -space) 
from  LASA,  was  overlaid  with  a  grid  containing  58  divisions  in  the  north- 
south  direction  and  48  divisions  in  the  east-west  direction.  The  edges  of 
the  grid  were  tangent  to  the  circle  corresponding  to  a  10-km/sec  velocity 
across  LASA.  Each  of  the  384  events  was  assigned  to  the  block  on  the  grid 
containing  its  epicenter.  The  average  traveltime  residual  (anomaly)  for 
each  subarray  was  computed  for  each  block  containing  one  or  more  events. 
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Several  adjacent  blocks  were  compared  using  the  mean  residual  computed  for 
each  subarray.  The  following  general  characteristics  were  observed  in  this 
traveltime  anomaly  analysis: 

•  Variations  in  traveltime  residuals  for  adjacent 
blocks  were  of  the  order  of  the  standard  devia¬ 
tion  of  the  data  in  each  block 

•  Variations  between  blocks  in  a  3-  to  5-block 
region  were  greatest  for  events  coming  from 
South  America 

•  Only  113  blocks  contained  events;  of  these, 
only  2  5  had  five  or  more  events 

•  Subarrays  on  the  E  and  F  rings  generally  had 
larger  means  and  standard  deviations  in  their 
residuals  than  did  the  inner  subarrays 


From  these  observations,  the  data  appeared  to  be  insufficient 
to  define  residual  correction  functions  covering  the  wavenumber  regions  of 
interest.  Until  sufficient  data  are  available  to  define  such  functions  with 
reasonable  confidence,  no  generalized  attempt  will  be  made  to  correct  the 
wavenumber  spectra  calculations  for  traveltime  anomalies. 


B.  MOVE  OUT  WINDOW  ANALYSIS 

A  simple  single-frequency  plane-wave  model  was  rsed  to 
analyze  the  spectral  window  due  to  a  "boxcar"  smoothing  function  applied 
to  a  crosspower  spectra.  When  a  crosspower  spectrum  is  formed  from 
two  finite  time  series,  some  type  of  smoothing  is  always  present.  The 
spectral  window  resulting  from  crosspower  smoothing  is  a  function  of  both 

the  time-series  length  (T)  and  the  signal  moveout  (T..,  the  signal  traveltime 

t  h  t  h  ^ 

from  the  i  to  the  j  sensor).  The  spectral  window  function  for  this  case  is 


H  (Af,  T  ) 

J 


sin  (n  Af  T  ij) 

(tTm  t..) 
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where  1  m  the  width  of  the  "boxcar"  smoothing  function.  The  smallest  if 
(least  smotrthing)  obtainable  with  a  time  aeries  of  length  T  is 

(  “  T 

This  results  in  the  following  relationship  between  data  length  T,  the  move  out 
between  channels  T.  ,  and  the  magnitude  of  the  crosspower  estimate 


where  ».  is  the  true  crosspower  between  channels  i  and  j.  The  crosspower 
ij 

estimate  will  then  have  a  magnitude  of  90  percent  of  the  true  crosspower 

magnitude  when  the  ratio  <  0.25. 

T 

Transform  gate  lengths  should  be  about  four  times  the  largest 
expected  moveout  across  the  array  if  crosspower  estimates  between  channels 
are  to  be  meaningful. 

The  importance  of  this  criterion  is  apparent  when  more  than 
one  signal  is  present.  In  the  case  of  two  plane  wavefronts  of  different  apparent 
velocity,  the  phase  of  the  crosspower  spcc’rum  har  ?n  ambiguous  interpreta> 
tion.  This  ambiguity  can  be  resolved  If  the  crosspower  magnitude  estimates 
between  pairs  of  sensors  are  close  to  the  true  crosspower  magnitudes. 

*♦ 

Currently,  the  effect  on  f-k  spectra  of  normalizing  the  cross* 
power  estimates  in  a  multiple-signal  environment  is  not  fully  understood; 
this  aspect  of  computing  high- resolution  f-k  spectra  is  the  object  of  further 
investigation. 
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C.  TRADEOFF  IN  FINITE -LENGTH  TRANSFORMS  OF  ARRAY  DATA 

A*  the  aperture  of  a  seismic  array  increases,  the  data  trans¬ 
form  gate  must  be  correspondingly  increased  to  keep  the  spectral  window 
from  deteriorating  the  information.  This  increase  in  transform  gale  reduces 
the  effective  signal-to-noise  ratio  when  the  signals  of  interest  are  transients 
within  the  gate.  This  decrease  in  signal-to-noise  ratio  can  be  offset  by 
adding  mure  sensors  as  the  aperture  is  increased,  assuming  that  signals 
of  the  same  minimum  apparent  velocity  are  applicable  to  both  arrays. 

Change  in  signal-to-noise  ratio  due  to  change  in  array  diameter  from  d 

k 

to  d^  and  change  in  number  of  sensors  from  Nk  to  N  can  be  expressed  as 


This  formula  indicates  that  doubling  the  number  of  sensors  and  douDling 
the  aperture  will  not  cnange  the  signal-to-noise  improvement  expected  in 
frequency  wavenumber  spectra  calculations.  At  LASA,  the  E  and  F  rings 
successively  double  the  aperture  of  the  previous  array  but  add  only  four  sen¬ 
sors  each.  The  available  signal-to-noise  improvement  is  then  decreased  when 
these  two  outer  rings  are  included  in  frequency-wavenumber  spectra  calculations. 

D.  SUMMARY 

F rom  this  section,  it  is  concluded  that  current  data  are  insuf¬ 
ficient  to  adequately  describe  the  amplitude,  traveltime,  and  waveform 
anomalies  at  LASA.  Additionally,  subarrays  on  the  E  and  F  rings  will  not 
be  used  in  frequency-wavenumber  spectra  evaluation.  These  subarrays  will 
be  used  for  teleseismic  event  detection  using  the  current  LASA  detection 
scheme. 


1-4 


solsnos  ssrvioss  division 


SECTION  II 


COMPARISON  OF  TRAVELTIME  RESIDUALS 

The  assumption  of  space  stationarity  necessary  for  computing 
meaningful  wavenumber  spectra  may  not  prove  to  be  a  valid  assumption  for 
large-diameter  arrays  such  as  LASA.  Departures  from  the  assumed  plane 
wavefront  of  constant  waveform  moving  at  constant  velocity  may  be  due  to 
two  primary  factors:  the  first  is  instrument  response  variations  and  should 
be  independent  of  wavenumber;  the  second  is  the  effect  introduced  by  different 
ray  paths  and  different  seismometer-to-earth  couplings.  Upper  mantle  in- 
hcmogenieties,  due  to  variations  in  thickness  and  composition,  will  probably 
be  a  function  of  wavenumber  (Appendix  B). 

At  least  part  of  the  second  factor  can  be  determined  empirically 
by  computing  the  residuals  caused  by  the  difference  between  theoretical  and 
actual  arrival  times  at  each  subarray  for  a  tele  seismic  P-wave  from  an  event 
of  known  epicenter.  These  computed  time  residuals  yield  the  phase  corrections 
used  when  computing  the  portion  of  a  wavenumber  spectra  corresponding  to 
that  epicentral  location.  In  this  report,  there  is  no  attempt  to  analyze  ampli¬ 
tude  anomalies  at  LASA. 

A.  WAVENUMBER  GRID  FOR  LASA 

Standard  wavenumber  spectra  computed  by  Texas  Instruments 
Incorporated  are  printed  out  on  an  alphabetically  coded  grid.  This  produces 
a  square  area  on  an  IBM-computer  printed  page.  Consistent  with  this  format, 
a  map  of  the  world  as  seen  tele seismicaliy  from  LASA  was  overlaid  with  a 
58-row  by  48-column  grid.  Rows  v/ere  aligned  in  the  east-west  direction. 

The  grid's  edges  were  tangent  to  the  circle  representing  the  of  epicenters 
which  would  theoretically  produce  a  10-km/sec  signal  velocity  across  LASA. 
Figure  II- 1  shows  the  wavenumber  grid  for  LASA. 
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B.  LASA  TRAVELTIME  ANOMAL/  COMPUTATIONS  FOR  BLOCKS  OF 
WAVENUMBER  SPACE 

The  wavenumber  plane  is  divided  into  a  grid  of  equal-sized 
blocks  to  correct  for  traveltime  anomalies  when  computing  frequency-wave¬ 
number  spectra  from  LASA  data.  Since  different  blocks  of  the  same  size  in 
wavenumber  space  cover  different  sizes  of  epicentral  area,  special  com¬ 
putation  is  needed  to  determine  the  anomalies  for  a  particular  block  from  the 
available  anomaly  information. 

The  TI  standard  wavenumber  printout  is  on  an  alphabetically 
coded  grid;  the  wavenumber  plane  is  divided  into  58  blocks  in  a  north-south 
direction  (k^)  and  into  48  blocks  in  the  east-west  direction  (k^).  The  edges 
of  the  grid  are  tangent  to  a  minimum- velocity  circle  of  10  km /sec.  A  par¬ 
ticular  block  in  the  grid  is  always  associated  with  the  same  area  on  the  sur¬ 
face  of  the  earth  as  seen  from  the  Montana  LASA. 

The  data  used  in  the  tabulation  include  3  84  events  with  date, 
latitude,  longitude,  and  arrival  time  (hours,  minutes,  and  seconds)  for  each 
subarray  on  punched  cards.  A  program  has  been  written  to  read  the  data  and 
perform  the  following  calculations. 

•  Epicentral  angle  (A)  and  azimuth  of  each 
event  for  each  subarray  are  calculated 
using  the  location  of  the  subarray  and  the 
following  formulas 

Event  epicentral  angle 

A  =  arc  cos  (Dl) 


where 

Dl  =  cos  (Bl)  •  cos  (Cl) 

+  sin  (Bl)  •  sin  (Cl)  •  cos  (Al) 
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and 

A1  =  (event  long.  )  +  (subarray  long.  ) 

B1  -  90.0°  -  (event  lat.  ) 

Cl  “  90.0°  -  (subarray  lat.  ) 

Azimuth  of  event 

Azimuth  =  Az  (event  in  NE  quadrant) 
Azimuth  =  180°  -  Az  (event  in  SE  quadrant) 
Azimuth  =  180°  +  Az  (event  in  SW  quadrant) 
Azimuth  =  360°  -  Az  (event  in  NW  quadrant) 


where 

|  arc  sin  (El)  | 


sin  (Al)  sin  (Bl) 
sin  (6) 

•  The  expected  traveltime  from  each  event  to 
each  subarray  is  determii  ed  using  a  second  - 
order  interpolation  between  the  1°  increments 
given  in  the  Jeffreys -Bullen  table 

•  The  traveltime  anomaly  of  each  event  for  each 
subarray  relative  to  subarray  A0  is  calculated 
using  the  formula 

AJ  *  (Tj  '  Hj)  ■  (Tao  -  hao) 


Az  = 

and 

El  = 


where 

T.  is  obser/ed  arrival  time  for  j**1  subarra* 
J 

H.  is  computed  traveltime  to  subarray 
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•  The  horizontal  velocity  of  each  event  at 
subarray  AO  is  interpolated  from  a  table 
using  epicentral  angle  & 

•  The  block  in  wavenumber  space  containing 
an  event  is  computed  using  the  horizontal 
velocity  and  azimuth  for  the  event 


After  the  preceding  five  steps  have  been 
completed  for  each  event,  the  wavenumber 
space  is  searched  block  by  block  and,  if  one 
or  more  events  are  found  in  a  block,  the 
traveltime  anomalies,  average  anomaly 
(AVERAGE),  and  standard  deviation  (SD^  are 
computed  for  each  subarray 
array, 


For  the  j1"  sub- 


AVERAGE  =  — 
N. 
J 


N. 

l 

£  A; 


1=1 


J 


SD 


i  N 

TTT  S  (Aj  - 


~i  1/2 


AVERAGE^ 


where 


th 

N.  is  number  of  events  in  j  block 
J 

A.  is  traveltime  residual  for  l  event  in 
J 

j1*1  block 


Table  II- 1  summarizes  the  distribution  of  the  384  usable  events 
over  the  113  blocks  within  which  their  epicenters  are  located.  The  calcula¬ 
tions  do  not  use  22  of  the  406  even  s  processed,  because  the  AO  arrival  time 
was  not  determined.  Figure  H-2  presents  distribution  of  the  events  over  the 
wavenumber  blocks. 
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Figure  II-2.  Number  of  EventB  in  Each  Block 
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Table  II- 1 

SUMMARY  OF  DISTRIBUTION  OF  USABLE  EVENTS 


No.  of  Events /Block 

No.  of  Blocks  with 
Given  No.  of  Events 

No.  of  Events 

1 

45 

45 

2 

20 

40 

3 

14 

42 

4 

9 

36 

5 

5 

25 

6 

5 

30 

7 

6 

42 

8 

2 

16 

9 

1 

9 

10 

1 

10 

12 

1 

12 

16 

2 

32 

17 

1 

17 

28 

1 

28 

Totals 

113 

384 

C.  COMPARISONS  OF  TRAVELTIME  RESIDUALS 

Comparison  of  the  traveltime  residuals  computed  for  the  tele- 
seismic  grid  (subsection  B)  was  made  using  blocks  in  the  Kazakh-Hindu  Kush, 
South  American,  and  Japan- Alaskan  regions  to  determine  the 

•  Difference  in  each  subarray  residual  for 
adjacent  blocks 

•  Difference  in  each  subarray  residual  at  2- 
and  3-block  distances 

•  Difference  over  a  3-block  region 

•  Difference  in  widely  separated  blocks  of 
the  same  general  region 

•  Comparison  of  block  residuals  to  previously 
computed  residuals  for  the  region  of  interest 

•  Portion  of  the  time  residual  due  to  the  dif¬ 
ferent  elevation  of  the  LASA  subarrays 
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These  comparisons  are  presented  in  Tables  II-2,  II-3,  and 
II-4.  The  Kazakh-Hindu  Kush  and  South  American  comparisons  were  made 
to  determine  the  differences  between  time  residuals  computed  for  adjacem 
blocks  in  wavenumber  space  and  previously  computed  time  residuals  for  these 
two  regions.  Two  sets  of  previously  computed  time  residuals  are  compared 
to  the  set  computed  for  this  report.  One  set  was  composed  of  time  residuals 
computed  by  TI  for  the  Large-Array  Signal  and  Noise  Analysis  Special  Report 
No.  1*  and  the  other  consisted  of  residuals  published  in  an  SDL  report. 

Many  of  the  events  used  in  the  SDL  residual  calculations  were  included  in  the 
data  set  used  in  the  calculations  for  this  report.  The  TI  Special  Scientific 
Report  No.  1  residuals  were  calculated  from  a  different  set  of  events. 

1.  Kazakh-Hindu  Kush  Comparison 


Only  one  event  in  the  data  card  library  had  an  epicenter  located 
in  block  17,23  -  the  Kazakh  region  of  the  USSR.  Residuals  computed  for  this 
event  compared  well  with  the  SDL-computed  residuals;  this  was  expected, 
since  the  SDL  Kazakh  event  ensemble  contained  this  event.  Residuals  com¬ 
puted  for  block  18,25  -  Hindu  Kush  region  -  were  identical  with  the  SDL- 
computed  residuals  due  to  a  nearly  identical  data  set.  Block  18,25  residuals 
agree  much  better  with  the  Kazakh  residuals  computed  for  TI  Special  Scien¬ 
tific  Report  No.  1  than  do  the  residuals  for  block  17,23;  this  indicates  that 
the  previous  Tl-calculated  residuals  contained  a  number  of  events  from  the 
Hindu  Kush  region. 


*  Texas  Instruments  Incorporated,  1967:  A  Study  of  the  Relative  Capability 
of  Large  and  Small  Seismic  Arrays  for  Event  Idennficat.on,  Large-Arr  y 
g  -  Rpt.  No.  1,  Contract  AF  33(657)- 1667b, 


Signal  and  Noise  Analysis,  Spec. 
20  Apr. 


“seismic  Data  Laboratory.  1966:  LASA  Traveltime  Anomalies  for  Various 
Epicentral  Regions,  ARPA  Order  No.  624,  *3  Sep. 
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KAZAKH  COMPARISON 
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SOUTH  AMERICAN  COMPARISON 


Tabic  II -4 


RELATIVE  ELEVATION  OF  SUBARRAYS  WITH  RESPECT  TO  AO 


Suba  r  ray 

Elevation^ 

(meters) 

Max.  Tirge.(msce 
at  4  km/ bcc) 

Diffcrercc  betwc 
17,  9  and  21.  ll 
(in  msec) 

AO 

— 

— 

-  -  ■ 

B  1 

+  10.0 

-  2.  5 

130 

B2 

-  50.  5 

12.  6 

140 

B3 

-  21.  9 

5.  5 

-  10 

B4 

-  27.  8 

7.0 

140 

Cl 

-  26.4 

6.  6 

150 

C2 

+  34.0 

-  8.  5 

190 

C3 

-  62.0 

15.  5 

70 

C4 

+  19.6 

-  4.9 

120 

D1 

+  14.  2 

-  3.5 

300 

D2 

-  83.  7 

20.  9 

80 

D3 

+  55.  1 

-13.  8 

190 

D4 

-  30.8 

7.  5 

130 

El 

-  58.  9 

14.  7 

180 

E2 

-134.  6 

33.6 

190 

E3 

+  16.  9 

-  4.2 

70 

E4 

+  58.  5 

-14.6 

10 

FI 

-  4.  3 

1.  1 

-  30 

F2 

+  9.  9 

-  2.5 

-160 

F3 

+  92.9 

-23.  2 

150 

F4 

-  37.0 

9.2 

£80 

+  is  above  AO 
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In  comparing  the  residual*  for  block  17,23  with  those  for 
block  18,25,  close  agreement  was  found  for  only  six  of  the  21  subarrays;  the 
remaining  15  subarrays  differed  by  at  least  one  LASA  sample  period  (50  msec). 
F4  was  the  only  subarray  in  the  E  and  F  rings  differing  by  less  than  50  msec, 
indicating  that  traveltime  residuals  for  Asia  changed  significantly  over  a  2- 
to  3- block  region.  Traveltime  residuals  for  the  E  and  F  rings  generally  have 
a  larger  mean  and  variance  than  the  residuals  computed  for  the  lower  rings. 

2.  South  American  Comparison 

Of  the  four  wavenumber  blocks  used  for  comparison  in.  the 
South  American  region,  three  were  adjacent  to  each  other  and  covered  the 
Northern  Argentina-Chile  border  region.  The  fourth  was  seven  blocks  away 
and  covered  part  of  the  Venezuelan  coast. 

Block  41,21  in  the  center  ot  the  3-block  group  showed  extremely 
close  correspondence  with  the  SDL  time  residual*  computed  for  this  region; 
this  was  expected,  since  most  of  the  events  used  were  common  to  both  com¬ 
putations.  Block  41,21  did  not  correspond  as  well  with  blocks  41,30  or  42,32 
as  it  did  with  the  SDL  Northern  Argentina-Chlle  residuals.  A  similar  change 
in  traveltime  residuals  for  adjacent  cpicentral  regions  in  South  America  can 
be  observed  in  the  SDL-computed  traveltime  residuals.  The  differences  in 
traveltime  residuals  within  the  3-biock  group  were  almost  as  great  as  the 
differences  between  biocks  41,31  and  44,37,  which  were  about  six  blocks 
apart.  This  indicates  degradation  in  wavenumber  resolution  when  several 
adjacent  blocks  within  the  South  America!.  region  are  corrected  with  the  same 
set  of  traveltime  residuals. 

3.  TravslUme  Residuals  Due  to  Subarray  Elevation 

As  shown  in  Figure  II-3,  the  nurber  of  subarrays  in  a  block 
needing  a  traveltime  correction  decreases  as  the  wavenumber  increases. 
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This  indicates  that  events  from  distant  epicenters  require  traveltime  cor¬ 
rections  for  a  greater  number  of  subarrays  than  do  events  of  closer  epicenter, 
even  though  the  distant  events  have  a  higher  apparent  velocity  across  the  array. 
It  was  thought  that  this  phenomenon  might  be  explained  by  a  difference  in 
elevation  between  each  subarray  and  the  reference  subarray. 

Table  II-4  summarizes  the  differences  in  elevation  and  the 
corresponding  traveltimes  in  a  medium  with  a  4-km/sec  propagation  velocity. 
These  computed  traveltimes  are  compared  with  the  difference  in  traveltime 
residuals  crmputed  for  blocks  17,  9  and  21, 18.  The  theoretical  traveltime 
anomalies  cue  to  elevation  are  insignificant  compared  with  anomaly  differences 
actually  experienced. 

D.  OBSER  VATIONS  IN  USING  ANOMALIES 

Experience  in  computing  high- re  solution  wavenumber  spectra 
at  LASA  indicates  that  regional  corrections  for  traveltime  anomalies  are 
usually  needed  to  detect  a  signal  using  the  full  LASA.  The  computed  subarray 
spectra  usually  can  detect  signals  without  anomaly  corrections.  Generally, 
the  effect  of  traveltime  anomalies  on  intermediate-size  arrays  (diametere 
from  20  to  60  km)  cannot  be  evaluated  using  the  currently  available  anomaly 
data.  These  effects  will  be  investigated  using  data  from  the  E3  extended 
subarray  and  from  a  LASA  configuration  consisting  of  the  A,  B,  C,  and 
D  rings  of  subarray  outputs. 

E.  CONCLUSIONS  ABOUT  ANOMALIES 

From  the  traveltime  anomalies  investigated,  the  following 
general  observations  were  made: 

•  Variations  in  traveltime  residuals  for 
adjacent  blocks  were  of  the  order  of  the 
standard  deviation  of  the  data  in  each 
block 
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•  Variations  between  blocks  in  a  3-  to  5-block 
region  were  greatest  for  events  coming  from 
South  America 

•  Differences  in  subarray  residuals  for  blocks 
separated  by  three  or  more  blocks  were 
greater  than  the  standard  deviation  of  the  data 
in  each  block 

•  Only  113  blocks  contained  events  and,  of  these, 
only  25  had  five  or  more  events 

•  Subarrays  in  the  LASA  E  and  F  rings  generally 
had  larger  means  and  standard  deviations  in 
their  calculated  residuals  than  did  the  inner- 
ring  subarrays 

•  Comparison  with  SDL  was  very  good  due  to 
common  data 

•  Relative  differences  in  subarray  elevation  had 
no  significant  effect  on  the  traveltime  residuals 

•  Comparison  with  previous  Tl-computed  residuals 
was  not  good  due  to  the  large  areas  used  in  pre¬ 
vious  calculations. 


Variability  of  the  mean  traveltime  residuals  was  small  enough 
to  permit  computation  of  average  residuals  for  multiple  block  regions.  This 

was  confirmed  in  Special  Scientific  Report  No.  1  where  a  large  region  was 

* 

considered  when  computing  residuals. 

The  data  currently  available  are  inadequate  for  the  design  of 
a  reasonable  scheme  to  correct  wavenumber  spectra  calculations  for  travel¬ 
time  anomalies.  This  was  apparent  from  the  wavenumber  distribution  of  the 
available  events.  Large  areas  of  interest  completely  lacked  events.  Before 
any  reasonable  correction  scheme  can  be  devised,  more  higher-quality  data 
must  be  available. 


* 

Op  cit 
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If  the  wavenumber  spectra  calculations  are  not  corrected  for 
traveltime  anomalies,  the  subarrays  on  the  LASA  E  and  F  rings  should  not 
be  included  in  the  calculations  due  to  their  generally  large  traveltime  resid¬ 
uals.  Other  arguments  for  not  including  these  subarrays  when  calculaving 
wavenumber  spectra  are  presented  in  the  following  sections. 
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SECTION  III 

SMOOTHING  MOVEOUT  ANALYSIS 


High-resolution  wavenumber  spectra  are  generated  by  using 
a  discrete  Fourier  transform  such  as  the  Cooley-Tukey  algorithm  to  trans¬ 
form  the  time -domain  data  to  the  frequency  domain  before  estimating  the 
auto-  and  crosspower  spectra.  This  section  uses  simple  models  to  analyze 
the  spectral  window  and  the  effective  signal -to-noise  improvement  experienced 
when  generating  high-resolution  spectra  from  LASA  data. 

By  estimating  the  crosspower  spectrum  of  two  sensors  in  a 
seismic  array  from  the  direct  Fourier  transform  of  the  individual  time  traces, 
a  spectral  window  effect  is  produced  in  the  estimate.  Under  the  simple  as¬ 
sumptions  of  uniform  plane-wave  propagation  across  the  array  and  a  boxcar 
smoothing  in  the  direct  transform,  this  spectral  window  is  expressed  as 


(3-1) 


where  T  is  the  length  of  the  data  trace  transformed  ana  t..  is  the  moveout  be¬ 
tween  the  i  and  channel.  If  the  energy  across  the  smoothing  window  is 
uniformly  distributed,  the  phase  of  the  crosspower  spectrum  is  unaffected  by 
the  spectral  window.  This  interpretation  (i.  e.  ,  processing)  of  a  crosspower 
spectrum  using  only  the  phase  information  becomes  ambiguous  when  there  is 
more  than  one  plane  wave. 

To  reduce  the  spectral  window  for  a  fixed  maximum  expected 
moveout  requires  that  the  transform  gate  length  be  increased.  Similarly,  to 
maintain  the  same  spectral  window  for  an  increased  maximum  moveout  also 
requires  an  increased  tran  form  gate  length.  This  increase  in  transform 
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gate  length  reduces  the  effective  signal-to-noise  ratio  for  the  transient  seis¬ 
mic  signal  being  considered.  In  both  cases,  the  decrease  in  effective  signal- 

to-noise  ratio  can  be  offset  by  adding  more  sensors  within  the  array. 


Although  the  following  development  treats  the  simple  cases 
just  discussed,  the  results  are  helpful  in  eliminating  gross  errors  when  gen¬ 
erating  crosspower  estimates  from  direct  transforms. 


A.  DEVELOPMENT  OF  WINDOW  FUNCTION 


Consider  two  sensors  S.  and  S.  separated  by  distance  d. .  in  a 

1  J  J 

seismic  array.  When  a  uniform  plane  wave  propagates  across  the  array  with 


an  apparent  velocity  of  v„  along  d„ ,  moveout  (traveltime)  t„  between  and 

S.  is 
J 


d.. 

t..  =  ju. 

1J  vij 


(3-2) 


We  let  6  be  the  angle  made  by  the  line  d. .  with  the  direction 
P  lJ 

of  propagation  of  the  plane  wave.  Then,  for  a  plane  wave  with  velocity  v  , 


apparent  velocity  v„  is 


(3-3) 


Moveout  t. .  is  then 
ij 


d. . 

t.  .  =  cos  0 

IJ  vp  P 


For  a  particular  v  ,  the  largest  moveout  which  can  be  ex¬ 


perienced  in  a  multisensor  array  occurs  between  sensors  with  the  widest 
separation  when  the  wavefronts  are  propagating  in  the  direction  of  separation 


(3-4) 


(i.e.  ,9  =0).  Later,  we  will  show  that  the  largest  expected  moveout  con- 

P 

strains  the  data  processing. 
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The  crosspower  spectrum  between  S,  and  S.  is  then 

i  J 


< p .  .  (£)  =  M  exp  (J  2ti  f  t,  ) 
ij  ij 


(3-5) 


Here,  we  assume  that  the  estimation  process  can  be  represented  as  a  boxcar 
smoothing  applied  over  an  interval 


f  s  f  s  f 
a  o  b 


(3-6) 


The  estimated  crosspower  at  frequency  f  is 


U 

0ij  =  /  eXP  <J  2"  £  ‘ij’  " 

a 


(3-7) 


Carrying  out  the  indicated  integration  gives 


0ij  =  2n  t. .  ™  -  (a)  [exP  <J  2"  ‘ij'b)  -  «P  U  2"  <3-8> 


Then,  the  real  and  imaginary  parts  of  Equation  3-8  are 


Re  0. .  =  — — 
rij  2tt  t. . 


M 


sin  2n  t. .  f,  -  sin  2tt  t 

(fK  -  f J  L  ij  b  ij 


ij  '  b  a' 


.0 


(3-9) 


and 


M 


Im  0.  .  =  - — - — — - TT  lcos  2tt  t. .  f  -  cos  2tt  t 

1J  2rT  j  (fb  "  fa>  L  iJ  a  ij 


(3-10) 
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By  employing  the  trigonometric  identities 


and 


sin  a  -  s 


cos  b  -  cos  a 


in  b  =  2  cos  (2y-^)  sin  (~ y^) 


*  2  8i"  (H*)  8in  (iTi) 


we  can  write  Equations  3-9  and  3-  10  as 


R'*ij 


M 


"*u  <fb  •  V 


and 


lm  = 


M 


ntij  <fb '  V 


jco.  [nt.,(fa+fb)]  .in  ["t4.  (fb  -  g]|  (3-H) 

j.in  [nt,.  (fb  +  fa)]  .in  [nt..  (£b  -  fa)]|  (3-12) 


The  magnitude  cf  the  crosspower  estimate  at  frequency  fQ  is  then 


I  ♦y  I  “  M 


sin  tt  t. .  (f, 
_LL  b 


TT  t,  .  (f 
ij  b 


and  the  phase  of  the  same  function  is 


(3-13) 


arg  0  =  tt  t  (£.+£)  (3-14) 

ij  ij  d  a 

In  general,  the  transform  algorithms  used  in  the  estimation 
process  have  a  constant-width  smoothing  function;  i.e.  ,  f.  -  f  =  Af  =  a  con- 

D  SL 

stant.  Thus,  Equation  3-13  can  be  written 


sin  tt  t. .  t>i 

M  - 4 — 

tt  t. .  Af 
ij 


(3-15) 
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B.  APPLICATION  TO  SEISMIC  ARRAYS 


The  smoothing-moveout  window  functions  are  present  in  the 
frequency-domain  processing  of  a  seismic  array  when  the  auto-crosspower 
matrix  is  being  estimated.  Usually,  the  processed  seismic  data  are  broad¬ 
band,  with  the  energy  nearly  uniformly  distributed  in  a  small  band  Af;.  in  this 
case,  the  crosspower  phase  estimated  by  Equation  3-14  will  be  a  very  good 
estimate  of  the  actual  phase,  but  the  smoothing-moveout  window  function  still 
will  be  present  in  the  magnitude  of  (he  crosspower  spectra  estimated  by 
Equation  3-15. 

The  effect  of  this  window  function  is  minin.ized  when  the  prod¬ 
uct  t.j  Af  is  minimized.  To  maintain  all  of  the  estimated  auto-  and  crosspower 
spectra  within  10  percent  of  each  other,  the  function 


sin  rrt  Af 

_ m 

TT  t  Af 
m 


a  0.  9 


(3-16) 


for  the  largest  expected  moveout  t  .  This  requires  that 

m 


t  Af  ^  - 
m  4 


(3-17) 


These  estimates  are  obtained  from  the  transforms  of  data  segments  T-sec 
long.  The  smallest  Af  obtainable  with  a  segment  of  length  T  is  approximately 


(3-18) 


This  implies  that 


T  ^  4  t 


The  length  of  the  data  segment  transformed  should  be  at  least 
four  times  the  largest  expected  moveout  across  the  array. 
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In  the  single  plane-wave  case  just  considered,  it  appears  that 
the  window  effect  could  be  circumvented  for  At  <  t,  .  <  tt  by  normalizing  the 
auto-  and  crosspower  spectra;  this  would  sacrifice  information  about  the 
waveform  amplitude  for  a  better  estimate  of  its  crosspower.  This  normal¬ 
izing  can  result  in  incorrect  interpretation  when  two  or  more  plane-wave  sig 
nals  of  different  velocities  are  present  in  the  transform  gate. 

For  two  plane  waves  with  moveouts  t^  and  t^,  respectively, 
the  crosspower  at  frequency  f  is 


(3-19) 


2  2 

0(f)  =  exp  j  2tt  ft  1 ,  +  M2  exp  j  2tt  ft2 
+  Mj  M2  exp  j  2tt  f  (tj  -  t2) 

+  M1  M2  exp  J  2tt  f  (^  ‘ 


where  and  M2  are  the  respective  amplitudes  of  the  waves.  The  estimated 
crosspower  can  be  expressed  by  integrating  each  term  of  Equation  3-19  over 
the  interval  f^  to  f^  and  operating  on  each  term  with  the  trigonometric  identi¬ 
ties  in  subsection  A.  The  estimated  crosspower  for  Af  =  f  -  f  is  then 

b  a 


sin  tt  t ^  Af 
_TTt  Af 


[1  exP  j  ntl  <£a  +  V 


(3-20) 


sin  tt  t2  Af 
."*2  Af 


exp  j  rr  t2  (f^  +  ffa) 


/sin  tt  (tj  -  t2)  Af^ 

I  _  Vi.  i.  \  A/ 

M.  M_  exp  | 

tt  (t  -  t  )(f  +  f  )] 

\  TT  -  *2I  01  /  *  *  l 

L  1  2  a  b  J 
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jl!  the  signals  are  wideband  and  the  energy  is  uniformly  dis¬ 
tributed  over  the  interval  Af*  th».n  (f  +  f,  )  is  very  close  to  2f  where  f  is 

a  b  o  u 

the  transform  frequency  point.  The  estimated  crosspower  then  contains 
every  term  of  the  true  crosspower,  with  each  term  multiplied  by  its  own 
window  function.  By  making  Af  small,  the  window  effects  can  be  reduced  and 
a  reasonably  good  estimate  obtained. 

C.  INCREASED  GATE  LENGTH’S  EFFECTS  ON  CROSSPOWER  ESTIMATES 

Current  applications  of  frequency-domain  processing  employ 
direct  transforms  of  the  time-domain  data  to  estimate  auto-  and  crosspower 
spectra.  Since  the  transform  algorithms  (such  as  Cooley- Tukey)  are  finite 
and  discrete,  their  outputs  are  functions  of  the  energy  rather  than  the  power 
contained  in  the  data  segment  transformed;  this  is  especially  true  in  seismic 
work  where  the  signals  being  processed  are  transients. 

Consider  the  transforms  X  and  X,  of  two  data  channels  con- 

a  b 

taining  signal  and  additive  noise.  Then, 

X  =  S  +  N 
a  a  a 


and 


xb  =  sb  +  Nb 


If  the  data  segment  is  of  length  T,  the  computed  transform 
points  will  have  a  frequency  resolution  limit  (spectral  window)  of  Afj  where 


M1  “  T 
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The  autopower  «pectrum  of  each  channel  is  then  estimated  as 
the  energy  density  spectrum 


—  *  r  *  *  v  *"} 

6  «  X  X.  -  S  S  +  SN  +  NS  +  NN 

a  ab  Laa  aa  aa  aaJ 


—  *  r  *  *  *  *1 

6.  *  X,  =  SL  Su  +  Su  Nl  +  Nu  Su  +  N 

rb  bb  Lbb  bb  bb  bbJ 

Similarly,  the  crosspower  spectrum  is  estimated  from  the  cross 
energy  density  spectrum 

—  *  r  *  *  *  *1 

6  L  =  X  XL  *  S  St  +  S  N.  +  N  S,  +  N  N 

ab  ab  Lab  ab  ab  abj 

In  these  expressions, 


sisi 


signal  energy  in  channel  i 


N.  N.  *  noise  energy  in  channel  i 

l  1 


5^  *  signal  energy  common  to  i  and  j  data  channels 


N,  N.  a  noise  energy  common  to  i  and  j  data  channels 

i  j 


For  transient  signals  with  an  effective  signal  duration  of  t  <  T, 

8 

the  band  Af  over  which  the  signal  energy  is  distributed  in  frequency  is  greater 
than  l/t^.  A  resolved  band  6f^  in  the  frequency  range  of  the  signal  energy 
will  contain  some  fraction  of  the  total  signal  energy  present  in  the  data  seg¬ 
ment.  If  the  data  segment  containing  this  signal  is  doubled  in  length  (2T),  the 
frequency  resolution  limit  df^  “0.5  &fj. 
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So,  if 


AfI  =  2Af2 


then 


Es  (Afj)  *  2  Eg  (Af2) 


Seismic  noise  is  generally  considered  to  be  the  output  of  a 
Gaussian  process.  Doubling  the  data  gate  length  will  then  double  the  total 
noise  energy  in  the  data  segment.  The  average  noise  energy  in  the  resolved 
band  A f2  using  a  double -length  (2T)  data  segment  will  then  be  the  same  as  the 
average  noise  energy  in  the  resolved  band  Af^  of  the  single-length  (T)  data 
segment;  i.  e.  , 

EN  (4fl>  ~  En  <W2> 

If  the  applicable  slgnal-to-noise  ratio  is  considered  to  be  the 
ratio  of  signal  energy  to  noise  energy  in  a  data  segment,  doubling  the  data  gate 
length  reduces  the  signal-to-noise  ratio  by  one-half: 

Eg  (A^2)  l  r.g  (Aft) 

V*f2>  =  2  [EN(Afl) 

This  is  true  whenever  signal  duration  t  is  less  than  the  data  gate  lengths 

s 

being  transformed. 

D.  LASA  PROCESSING  CONSIDERATIONS 

Results  of  the  previous  analysis  on  the  computation  of  high- 
resolution  wavenumber  spectra  for  the  Montana  LASA  may  be  applied  by  con¬ 
sidering  the  beamsteered  improvement  in  signal-to-npise  power.  This  ap¬ 
proach  is  useful  because  the  computation  of  standard  wavenumber  spectra  is 


IH-IO 


science  services  division 


the  frequency -domain  equivalent  to  time -shift-and-sum  beam-forming  fol¬ 
lowed  by  square-law  detecting.  In  the  frequency-domain  processing  con¬ 
sidered,  all  time -domain  data  are  transformed  to  the  frequency  domain  be¬ 
fore  the  performance  of  any  beamsteering  or  other  processing.  The  relative 
signal-to-noise  gain  for  various  LASA  configurations  will  then  give  a  general 
picture  of  the  manner  in  which  configuration  and  moveout  constrain  the  pro- 
ces sing. 

To  limit  the  extent  of  this  comparison,  the  array  configurations 
include  all  of  the  subarrays  within  and  on  the  LASA  ring  under  consideration. 
Thus,  there  are  five  array  configurations  for  the  B  through  the  F  rings.  The 
fignals  of  interest  have  greater  than  10-km/sec  apparent  velocity  across  the 
array,  so  the  maximum  expected  moveout  for  each  configuration  is 

d. 

t  =  _ J _ 

m.  10  km/sec 
J 

where  d.  is  the  diameter  of  the  j1-*1  ring. 

J 

If  the  array  eu  irged  by  adding  all  of  the  subarrays  on  the 

next  ring  out,  the  transform  gate  length  T.  mu«t  be  increased  to  maintain 

J 

the  same  spectral  window;  this  increases  the  transform  gate  length  to 


The  loss  in  effective  signal-to-noi se  ratio  due  to  the  increase  in  gate  length 
is  then 


T.  d. 
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As  the  array  is  expanded,  more  sensors  are  added.  Assum¬ 
ing  that  the  seismic  noise  field  is  uncorrelated  from  subarray  to  subarray, 
the  beamsteered  signal-to-noise  power  gain  of  the  j**1  array  over  a  single 
subarray  is 

rJ’j 

G.  =N,  =  number  of  subarrays  in  jv  configuration 
J  J 

The  total  gain  (or  loss)  obtained  by  adding  one  ring  to  the  j1*1  array  configura 
tion  is  then 


v.j  ■  v*  U~ 

\  J  v 


d,  \  /N. 


,VAnj 


As  successively  larger  rings  are  added  to  the  beamsteering 
process,  the  signal-to-noise  improvement  predicted  for  LASA  increases  with 
the  increase  in  the  number  of  subarrays  and  decreases  with  the  increase  in 
aperture.  This  function  is  plotted  in  Figure  III-2. 

This  analysis  indicates  a  serious  drawback  in  using  rings  out¬ 
side  the  C  ring  for  frequency-domain  sijnal  processing.  To  include  the  D 
ring  reduces  the  beamsteer  gain  20  percent,  while  the  E  and  F  rings  reduce 
the  beamsteer  gain  50  and  60  percent,  respectively.  These  decreases  result 
from  doubling  the  array  aperture  while  not  doubling  the  number  of  subarrays. 

In  the  light  of  this  analysis,  the  use  of  subarrays  on  the  E  and 
F  rings  when  computing  high- resolution  wavenumber  spectra  for  LASA  does 
not  appear  advantageous. 
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SECTION  IV 
CONCLUSIONS 

From  this  study,  there  are  two  major  conclusions: 

•  Current  data  are  insufficient  to  define  a 
scheme  adequately  to  correct  wavenumber 
spectra  calculations  for  traveltime  anomalies 

•  Subarrays  on  the  E  and  F  rings  of  LASA  will 
not  be  included  in  high- re  solution  f-k  spectra 
calculations 

Subarrays  on  the  E  and  F  rings  of  LASA  generally  exhibit 
larger  traveltime  residuals  and  less  waveform  similarity  than  do  subarrays 
of  the  inner  rings.  Current  data  appear  to  be  inadequate  to  describe  these 
anomalies  with  any  degree  of  certainty.  With  this  lack  of  knowledge  of  the 
anomaly  mechanism,  coupled  with  the  spectral  window  considerations,  the 
use  of  subarrays  on  the  E  and  F  rings  of  LASA  is  unadvisable. 


Texas  Instruments  Incorporated,  1967:  Short-Period  Signal  Waveform  at 
LASA,  Large-Array  Signal  and  Noise  Analysis,  Spec.  Rpt.  No.  8,  Con¬ 
tract  AF  33(657)- 16678,  1  Aug. 
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APPENDIX  B 


CORRECTIONS  FOR  DEPARTURES  FROM  SPACE  STATIONARITY 

WHEN  COMPUT'NG  WAVENUMBER  SPECTRA 

The  assumption  Ox  space  stationarity  necessary  for  the  com¬ 
putation  of  meaningful  wavenumber  spectra  may  not  prove  to  be  valid  for 
large-diameter  arrays  such  as  LASA.  Departures  from  the  assumed  plane 
wavefront  of  constant  waveform  moving  at  constant  velocity  may  be  due  to 
two  primary  factors:  the  first  is  instrument  response  variation.,  and  should 
be  independent  of  wavenumber;  the  second  is  the  effect  introduced  by  different 
crustal  paths  and  different  seismometer -to-earth  couplings.  Upper  mantle  in- 
homogenieties,  due  to  variations  inthickness  and  composition,  will  probably  be  a 
function  of  wavenumber.  If  the  total  effect  of  these  two  factors  can  be  de¬ 
termined  theoretically  or  empirically,  corrections  for  departures  from  space 
stationarity  may  be  easily  made  as  follows. 

Let  H  (w,  k)  be  the  transfer  function  of  the  filter  which  equal- 
.  |  J 

izes  the  j  seismometer  (or  subarray)  to  the  r*"^1  or  reference  seismometer 
(or  subarray).  To  compute  the  power  density  at  frequency  w  and  wavenumber 
k,  the  matrix  $  of  auto-  and  crosspower  spectra,  should  be  premultiplied  by 
T  and  postmultiplied  by  T*  as  shown: 


B-l 
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m 


When  computing  high-resolution  wavenumber  spectra  (esti¬ 
mating  the  signal  at  location  j)  for  the  spatially  random  signal  model  case, 
the  matrix  equation  solved  for  the  filter  responses  F^w),  F2(w), . . 


F  (w)  is 
n 


where 


»  (w)t_ia=L 

"  |H,(w)|2 


»21(w) 


$  (w) 

N1  ' 


K(w) 


|  H  (w)  | 
J 


**  =  L  ♦ 


f12(w) 


$1N(W) 


f  (w)f  K(W^ 
|H2(w)  | 


Ww) 


K(w) 


NN  |hn(w)|2 


F  *(w) 


f2V) 


F*  - 


F*(w, 
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n 


To  correct  for  departures  from  space  stationarity,  we  should 
solve  the  matrix  equation 


T  £  T*  F  *  =  H*(w)  T f  (B-l) 

c  J 

or 

F  *  =  H.*(w)  [T*]“ 1  £  T_1T*  (B_2) 

c  J 

F *  =  H*(w>  [T*]"1  £  1  ft  =  H.*(w)  [T*]_1F* 

c  .1  J 

Thus,  the  corrected  filter  responses  may  be  obtained  from  the 
uncorrected  filter  responses  from  Equation  B-2  more  simply  than  through 
the  solution  of  Equation  B-l. 

A  similar  result  is  obtained  for  the  multichannel  Markov  case 
for  the  correction  to  the  spatial  prediction  filter  responses.  The  equation 
to  be  solved  for  the  filter  to  predict  the  N+lth  channel  now  is 

T$T*F  *  -  H  *(w)Tr 

c  N+ 1 

wheie 


Therefore, 


F .*  =  [t*]'1  »'1t'1tt  *(w>  =  HM,  ♦(w)LT*]"1  t1  r 


Fc*  =  Vi‘(w)[TT1f* 


where 


f*  =  $_1r 
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